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In the last decade, computational studies of liquid water have mostly concentrated on ground 
state properties. However recent spectroscopic measurements have been used to infer the structure 
of water, and the interpretation of optical and x-ray spectra requires accurate theoretical models 
of excited electronic states, not only of the ground state. To this end, we investigate the electronic 
properties of water at ambient conditions using ab initio density functional theory within the gener- 
alized gradient approximation (DFT/GGA), focussing on the unoccupied subspace of Kohn-Sham 
eigenstates. We generate long (250 ps) classical trajectories for large supercells, up to 256 molecules, 
from which uncorrelated configurations of water molecules are extracted for use in DFT/GGA cal- 
culations of the electronic structure. We find that the density of occupied states of this molecular 
liquid is well described with 32 molecule supercells using a single k-point (k = 0) to approximate 
integration over the first Brillouin zone. However, the description of the density of unoccupied states 
(u-EDOS) is sensitive to finite size effects. Small, 32 molecule supercell calculations, using the T- 
point approximation, yield a spuriously isolated state above the Fermi level. Nevertheless, the more 
accurate u-EDOS of large, 256 molecule supercells may be reproduced using smaller supercells and 
increased k-point sampling. This indicates that the electronic structure of molecular liquids like 
water is relatively insensitive to the long-range disorder in the molecular structure. These results 
have important implications for efficiently increasing the accuracy of spectral calculations for water 
and other molecular liquids. 

PACS numbers: 71.15.Dx,71.15.Mb,78.40.Dw,78.40.Pg 



I. INTRODUCTION 

Our understanding of the structure of water in its 
many phases is fundamental to research in fields as 
diverse as biochemistry, cellular biology, atmospheric 
chemistry, and planetary physics. The standard exper- 
imental approaches to the analysis of crystal structures 
- x-rajii and neutron^ diffraction - can provide detailed 
information on ordered phases of ice, but only indirect, 
and often limited, information on amorphous ice and liq- 
uid water. Nonetheless, the results of these experiments 
have permitted major advances in the understanding of 
the properties of water in the last 30 years; in addition 
they made possible the development of simple classical 
potentials to describe water in the condensed phase, thus 
enabling molecular dynamics simulations to probe its dy- 
namical properties. 

Recently, more sophisticated ab initio electronic struc- 
ture approaches have become available to perform simu- 
lations of water & Ab initio molecular dynamics based on 
density functional theory 4,5 (DFT) has been extensively 
used to study water and solvation processes. It is not yet 
fully understood how accurate DFT is&£ - and, in partic- 
ular, how accurate the various gradient corrected func- 
tionals&Siifl are - in describing the structural and diffu- 
sive properties of liquid water. Furthermore, the quanti- 
tative influence of the inclusion of proton quantum effects 
in ab initio simulations remains to be established^ How- 
ever, recent results reported in the literature have shown 
that qualitative features of hydrogen bonding in liquid 
water are correctly accounted for by ab initio molecular 
dynamics (MD) based on DFTi 6 : 11 : 12 

The majority of ab initio studies, to date, have concen- 



trated on structural properties of liquid water, and, only 
recently, have first-principles calculations of absorption 
spectra been carried outiiii^ii 5 ^* 1 ^ 1 ^* 1 ^ These compu- 
tations were motivated by new spectroscopic results ob- 
tained using x-ray Raman spectroscopjiiS*22iSi and x-ray 
absorption spectroscopy 15,17 ' 18,2 ^ for water and ice; and 
x-ray absorption spectroscopy for liquid water jets^*^ 
The experimental spectra - providing direct information 
on the electronic transitions from atomic core levels to 
excited states - have been used to infer information on 
the structural properties of the fluid. In particular, the 
number of hydrogen bonds and the details of the hydro- 
gen bonded network in the fluid have been inferred from 
experiment and compared to those of ice. A standard 
procedure has been to use DFT in the gradient corrected 
approximation to compute absorption spectra of selected 
snapshots representing the liquid - within some approx- 
imation for the description of the core-hole interaction. 
Those snapshots producing spectra in agreement with ex- 
periment have then been considered to be representative 
of the "correct" water structure, as probed experimen- 
tally. While this may well be a viable and straightfor- 
ward approach in the absence of major approximations 
in the theory used to compute spectra, it becomes a much 
more complex interpretative tool in the presence of ap- 
proximations in the theory, in particular approximations 
regarding the core-hole interaction with the excited elec- 
tron. 

In order to provide a clear interpretation of measured 
and computed spectra, an important prerequisite is to 
fully understand the electronic structure of the fluid, as 
described within DFT. The purpose of this article is to 
provide a detailed description of the electronic structure 
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of liquid water using Kohn-Sham density functional the- 
ory, and in particular to understand in detail the unoccu- 
pied subspace of Kohn-Sham eigenstates - especially the 
conduction band minimum, which will have a great in- 
fluence in determining the onset of computed absorption 
spectra. 

We shall show that the electronic density of states 
(EDOS) of liquid water, near the conduction band min- 
imum, is particularly sensitive to finite size effects. We 
find that for small, 32 molecule supercell calculations of 
the DFT electronic structure of liquid water, approximat- 
ing Brillouin zone integration by sampling the L-point 
only, leads to an EDOS where the lowest unoccupied 
state is separated from the rest of the conduction band 
by ~ f .5 eV. This is consistent with previous DFT cal- 
culations of the electronic structure of water. 

In order to understand this peculiar feature of the 
EDOS of water, we first examine the isolated dimer 
and the ordered hexagonal phase of ice Ih. The Kohn- 
Sham eigenstate energies of the dimer indicate a separa- 
tion in energy of the lowest unoccupied molecular orbital 
(LUMO) from the rest of the unoccupied states. This is 
apparently consistent with the separation of such states 
in liquid water. However, water possesses a hydrogen 
bonded network more similar to ice than to the isolated 
dimer. An examination of the band structure of the ice Ih 
crystal indicates large dispersion in the conduction band 
states. In particular, there is a large separation (~ 3 eV) 
between the two lowest conduction bands at the T-point. 

Inspired by these observations, we examine the conver- 
gence of the EDOS of liquid water with respect to k-point 
sampling. We compute the EDOS using several uncorre- 
cted configurations of water molecules generated using 
molecular dynamics at ambient conditions with a clas- 
sical potential. Examining the EDOS of several liquid 
water configurations, we find that the occupied portion 
is converged using just one k-point. This implies that the 
results of previous DFT studies of water in the electronic 
ground state would be unaffected by increased k-point 
sampling. However, for the unoccupied portion of the 
EDOS, the separation between the two lowest conduc- 
tion bands is merely an artifact of poor k-point sampling 
of the Brillouin zone. We test whether the EDOS com- 
puted using a 32 molecule supercell with 8 k-points is an 
accurate approximation of the EDOS computed using a 
256 molecule supercell with 1 k-point. Such calculations 
would yield identical results for a periodic crystal; how- 
ever, for a disordered phase like liquid water this equiv- 
alence is destroyed. Yet, at 300 K the differences are 
small, and we find that increased k-point sampling in 
smaller supercells is indeed an accurate approximation 
of the EDOS of much larger liquid water systems. This 
indicates that our converged, k-point sampled EDOS of 
liquid water is accurate and that the presence of an iso- 
lated LUMO of water is spurious. 

Our analysis has important consequences for the cal- 
culation of spectral properties (both optical and x-ray) of 
pure water and materials in aqueous solvation, since the 
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FIG. 1: Occupied and unoccupied Kohn-Sham eigenenergies 
in the neighborhood of the Fermi energy (HOMO), computed 
using DFT/PBE for sequential molecular configurations of 
32 water molecules taken from a 20 ps Car-Parrinello MD 
simulation at 300 K, using only the F-point. There exists a 
clear separation in energy between the LUMO and the rest of 
the conduction band. 



accuracy of these quantities is dependent on a reliable 
description of the unoccupied EDOS. 



The outline of our work is as follows: Before present- 
ing our results, we discuss some possible misunderstand- 
ing of the conduction band in water, which exists in the 
literature (Section [H]). To this end, we build some in- 
tuition on the electronic structure of liquid water based 
on an isolated water dimer and a crystalline phase of 
ice (Sections lllll and II VI respectively). We efficiently 
generate a large number of molecular configurations of 
liquid water, for use in DFT calculations, using classi- 
cal molecular dynamics, as indicated in Section In 
Sections I V II and I V If I we make use of these MD trajecto- 
ries to illustrate that, while the L-point approximation is 
valid for Brillouin zone integration in the occupied sub- 
space of Kohn-Sham eigenstates in liquid water, this is in 
fact a very poor approximation for the unoccupied states. 
To indicate the accuracy of the electronic structure com- 
puted with increased k-point sampling, we compare in 
Section IVIIII with the EDOS of larger supercell calcu- 
lations. We outline in Section ITXl the implications that 
such improvements in the description of the electronic 
structure will have on calculations which make use of the 
unoccupied EDOS, and in Section [X] we show that our 
approach is quite efficient as a means of increasing the 
accuracy of these calculations. Finally, we give our con- 
clusions in Section |XH 
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II. PREVIOUS DESCRIPTIONS OF THE 
LOWER CONDUCTION BAND OF LIQUID 
WATER 

The earliest ab initio calculations on liquid water- 
reported the existence of a delocalized lowest unoccu- 
pied molecular orbital (LUMO). This LUMO comprised 
some molecular a* character on the oxygens of each water 
molecule and a delocalized tail which extended through- 
out the system in the intcrmolecular volume, avoiding 
hydrogen-bonds. In addition, the reported electronic 
density of states (EDOS) for a supercell of 32 heavy- water 
molecules (D2O), calculated within DFT using the Becke- 
Perdew gradient corrected exchange correlation func- 
tional£i££^£ reveals an unexpected, isolated state at the 
bottom of the conduction band. This is perhaps coun- 
terintuitive, as one might expect that a molecular system 
would possess bands of almost degenerate states arising 
from weak coupling of the original isolated molecular or- 
bitals; therefore, the conduction band minimum would 
be the lowest in energy of a band of states produced by 
coupling of the LUMOs of single water molecules. 

Later workSS displayed the Kohn-Sham eigenenergies 
of water, as a function of time during a Car-Parrincllo 
molecular dynamics 3 simulation, at densities lower than 
the ambient density. This showed the isolated peak above 
the Fermi energy of liquid water to be caused by a single 
Kohn-Sham eigenstate, whose separation from the rest 
of the conduction band grows with water density. The 
authors of this work infer that this "delocalized LUMO 
can be regarded as the precursor of the solvated elec- 
tron." At these low densities they state that this LUMO 
has more delocalized character than those other virtual 
orbitals higher in energy, which possess more molecu- 
lar character. Subsequent workSi simulated the occupa- 
tion of this LUMO with an excess electron, performing 
Car-Parrinello MD at ambient conditions, leading to the 
formation of a localized region of charge at intervals of 
approximately 0.1 ps mediated by the breaking of hydro- 
gen bonds. The calculated optical absorption from these 
simulations showed good agreement with hydrated excess 
electon experiments— 

Since then, a large body of work on water and solvated 
systems has generated similar time-dependent EDOS, 
all exhibiting a clear separation of the LUMO of wa- 
ter from the rest of the conduction band. ^i i 32 4 33 i 34 
For clarity, we report in Fig. ^ a similar set of Kohn- 
Sham eigenenergies as a function of time for a Car- 
Parrinello MD simulation of 32 water molecules at am- 
bient density (p=1.0g/cm 3 ) and temperature (300 K). 
The Car-Parrinello trajectory was generated using the GP 
code 35 and the Kohn-Sham eigenenergies of the unoccu- 
pied states for the configurations taken from this trajec- 
tory were computed using ABINIT. 36 Both sets of calcu- 
lations used DFT in the generalized gradient approxima- 
tion (GGA) of Perdew, Burke and Ernzerhof (PBE)i2I 
These were planewave, norm-conserving pseudopoten- 
tialf^ electronic structure calculations performed in a 
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FIG. 2: The water dimer in its ground state geometry com- 
puted using DFT/PBE. Left: Isosurfaces of the probability 
densities of the first two unoccupied Kohn-Sham eigenstates 
above the Fermi energy - LUMO and LUMO+1. Water 
molecules are indicated using red (oxygen) and white (hy- 
drogen) ball-and-stick models. Right: A comparison of the 
Kohn-Sham eigenenergies of the monomer and the dimer indi- 
cating a separation in energy between LUMO and LUMO+1. 
Note that the energies in both cases have been shifted to align 
the HOMO with zero. 



supercell under periodic boundary conditions. We trun- 
cate the planewave basis using a kinetic energy cut-off 
of 70 Ry and the oxygen pseudopotential is non-local in 
the s angular momentum channel. The Brillouin zone 
integration is approximated using a single k-point, k = 
(the T-point). This is a standard approximation for rea- 
sonably large supercells of insulating systems and is the 
same approach used in previous work, which led to EDOS 
of water similar to that of Fig. with an isolated LUMO 
state. 

The presence of such an isolated LUMO is not appar- 
ent from the experimentally measured optical spectrum 
of water (2242i which displays no noticeable peak near the 
onset of absorption. This, however, does not discount the 
presence of such a peak in the unoccupied EDOS, which 
could be diminished in the spectrum by reduced oscilla- 
tor strengths associated with optical transitions from the 
top of the valence band. In the following sections we shall 
examine the electronic structure of the water dimer (the 
smallest hydrogen bonded water system) and of hexago- 
nal ice (a fully saturated hydrogen bonded phase). From 
this analysis we shall develop an improved picture of the 
unoccupied EDOS of liquid water and show, ultimately, 
that there is no isolated peak in the EDOS at the bottom 
of the conduction band. 



III. THE WATER DIMER 

Past justifications for the presence of an isolated 
LUMO in the u-EDOS of liquid water emphasized com- 
parisons with the molecular orbitals of the water dimer. 
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We used DFT/PBE to compute the electronic structure 
of the water monomer and dimer. Note that all subse- 
quent DFT calculations outlined in this work were per- 
formed using the PWSCF package^ and with a planewave 
cut-off of 85 Ry, unless stated otherwise . Two Kohn- 
Sham eigenstates of the dimer, the LUMO and the Kohn- 
Sham eigenstate just above it in energy (LUMO+1), are 
indicated by isosurfaces of their probability densities in 
Fig. [21 It is observed by comparison of the Kohn-Sham 
eigenenergies of the dimer and monomer that the LUMO 
is lower in energy than the LUMO+1 as a consequence 
of being localized on the dangling hydrogen bonds in this 
system. Similarly, in the case of liquid water, the LUMO 
is seerj^ to possess a similar a* component on the oxy- 
gens in the system and also a significant probability den- 
sity in those regions where the hydrogen-bond network is 
disrupted. 



IV. BAND STRUCTURE OF ICE Ih 

Liquid water is perhaps more akin to ice than a gas 
phase dimer. In ordered phases of ice, all hydrogen 
bonds are saturated (i.e., four hydrogen bonds per water 
molecule) and, molecular dynamics simulations of liquid 
water at ambient conditions indicate that the average 
number of hydrogen bonds per molecule is between 3 
and 4, depending on the particular definition of a hydro- 
gen bond4£ Therefore, before analyzing liquid water in 
detail, we choose to analyze the electronic structure of 
ice. 

Previous work^iM on the band structure of cubic ice 
(Ic) and the density of states of hexagonal ice^ (Ih) 
gives no indication of peculiarity at the conduction band 
minimum. However, the band structures indicate sig- 
nificant dispersion in the bands just above the Fermi 
energy. Using the hexagonal unit cell of ice Ih as pro- 
posed by Bernal and Fowled, we use DFT/PBE cal- 
culations to generate the structural parameters within 
the Born-Oppenheimer approximation at a pressure and 
temperature of ~ 0.1 MPa and K, respectively. Pro- 
ton disorder is not considered while we use the primitive 
unit cell containing 12 water molecules with all hydrogen 
bonds passivated (Fig. |3J. This structure has a density 
of 1.002 g/cm 3 , corresponding to hexagonal lattice pa- 
rameters of a — 7.59 A and c = 7.18 A. The use of 
the PBE functional is justified by previous work on ice 
at various pressures^ We present the band structure in 
Fig. 2] The self-consistent electronic charge density is 
computed using 8 k-points in the Brillouin zone and all 
eigenenergies in the band structures computed non-self- 
consistently using this charge density and its associated 
Kohn-Sham potential. Our band structure is similar to 
that reported by Hahn et al. ^ where they employ a cubic 
cell with some proton disorder, and we use the primitive 
hexagonal cell with no proton disorder. 

We note immediately that the large degree of disper- 
sion in the unoccupied subspace leads to a separation of 




FIG. 3: The ice Ih structure with isosurfaces of the probability 
density (green) of the Kohn-Sham eigenstate just above the 
Fermi energy at the F-point: (a) the oxygen a* component; 
and (b) the delocalized, hydrogen-bond-avoiding component, 
corresponding, respectively, to ~ 5% and ~ 30% of the in- 
tegrated density of the state. Water molecules are indicated 
using red (oxygen) and white (hydrogen) ball-and-stick mod- 
els. Hydrogen bonds indicated as dashed lines. 



the lowest conduction bands similar to that between the 
LUMO and LUMO+1 in previous calculations for water. 
This separation is particularly large at the T-point. It 
is clear that considering only the band structure at the 
T-point in Fig. 0] would lead to a similar EDOS as that 
outlined for liquid water in Fig.^ Also, the Kohn-Sham 
eigenstate at the bottom of the conduction band is de- 
localized, with a* character on the oxygen atoms and 
an avoidance of hydrogen bonds, similar to the LUMO 
in previous water calculations. Furthermore, the sepa- 
ration between the lowest two conduction bands at the 
T-point is 3.0 eV. This large energy is consistent with the 
trends indicated by Boero et al. for the variation of this 
separation with density.^ The discrepancy between this 
large separation in ice Ih and that of ~ 1.5 eV in ambient 
liquid water (which has a similar density) is likely due to 
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FIG. 4: The band structure of ice Ih computed using 
DFT/PBE for the 12 molecule unit cell, with structural pa- 
rameters optimized to reduce the pressure to ~ 0.1 MPa. The 
electronic charge density is determined using 8 k-points in the 
first Brillouin zone, (see Fig. [3] and text). 

the saturation of all hydrogen bonds in the former, lead- 
ing to a more delocalized state and consequently more 
dispersion at the T-point. 
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FIG. 5: The oxygen-oxygen radial distribution function of 
liquid water at ambient conditions determined experimen- 
tally- (red) and from a molecular dynamics simultion using 
the TIP4P classical potential. 

dial distribution function of water at ambient conditions. 

We used the gromacs molecular dynamics package 53 " 52 
to generate long trajectories for large supercells of water 
molecules in the NVT ensemble. We choose the density 
to be 0.996 g/cm 3 and the temperature to be 300 K. MD 
simulations for all supercell sizes reported here are begun 
using a box of water molecules cut from a large (2048 
molecule) equilibrated sample. The molecules in this su- 
percell are provided with a Boltzmann velocity distribu- 
tion consistent with a temperature of 300 K and allowed 
to equilibrate for 50 ps. A further 200 ps is evolved for 
sampling purposes. All configurations used for DFT cal- 
culations are separated by 20 ps. Previous work indicates 
that this separation is of the same order of magnitude as 
the structural correlation time. «£ 



V. CLASSICAL TIP4P TRAJECTORIES FOR 
LIQUID WATER 

In an attempt to facilitate the efficient reproduction 
of our results, we use configurations of water molecules 
generated using a classical potential. This removes the 
computational expense of generating such configurations 
using ab initio MD, and also provides configurations 
with structures in closer agreement with experiment, at 
least for such measures as the radial distribution func- 
tions and diffusion coefficient. Recent, careful, ab initio 
DFT/GGA MD simulations have been shown to produce 
more structured radial distribution functions in compar- 
ison with experiment&Lii We use the TIP4P four-site 
models^, which has been shown recently to have a wide 
range of transferability for water in various condensed 
phases^ We see in Fig. [5] that TIP4P approximates 
well the experimentally determined^ oxygen-oxygen ra- 



VI. "BAND STRUCTURE" OF LIQUID WATER 

We provide in Fig. the DFT/PBE calculated "band 
structure" for a representative cubic supercell of 32 wa- 
ter molecules extracted from the equilibrated section of 
a TIP4P classical trajectory. We recognize that the con- 
cept of a band structure has no meaning for an aperi- 
odic system, however the band structures of the periodic 
approximations to the true disordered systems provide 
information about the electronic structure of the liquid. 
Convergence of the electronic charge density with respect 
to k-point sampling in self-consistent calculations for this 
system indicated that using just the T-point is a valid ap- 
proximation rather than attempting a full integration of 
the Brillouin zone. This fact is clearly evidenced by the 
lack of significant dispersion in the occupied subspace of 
this liquid water configuration. We computed the Kohn- 
Sham eigenvalues at each k-point in this band structure 
using the effective potential derived from the electronic 
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FIG. 6: Computed DFT/PBE "band structure" for one rep- 
resentative 32 molecule supercell of liquid water taken from a 
TIP4P trajectory. The electronic charge density is converged 
using the T-point only. 



charge density computed using the T-point approxima- 
tion. 

Examination of the unoccupied bands in this small liq- 
uid water supercell reveals large dispersion, particularly 
at the F-point. This indicates that the u-EDOS of this 
periodic system will be sensitive to the amount of k-point 
sampling employed. Clearly, using the EDOS at the T- 
point will lead to the familiar separation of states at the 
bottom of the conduction band as seen in Fig. ^ How- 
ever, we can see now that increasing k-point sampling 
will add to the EDOS in the region that has traditionally 
separated the LUMO and LUMO+1 in T-point calcula- 
tions of liquid water using supercells containing 32 water 
molecules. We examine this convergence in the next sec- 
tion. 



VII. CONVERGENCE OF THE EDOS WITH 
K-POINT SAMPLING 

We take 10 uncorrelated molecular configurations of 
32 water molecules of liquid water, at ambient condi- 
tions, at intervals of 20 ps from a TIP4P trajectory. For 
each of these we perform DFT/PBE electronic structure 



calculations, increasing the k-point sampling of the first 
Brillouin zone. The indicated number of k-points refers 
to the number of points in a uniform grid about k = 0. 
For example, 8 k-points implies a 2 x 2 x 2 grid. How- 
ever, using the symmetry of our cubic supercells we re- 
duce the actual number of k-points used in the calcu- 
lation, reweighting each appropriately in the sum which 
approximates a complete Brillouin zone integration. In 
each of the configurations examined we find that 64 k- 
points is sufficient to converge the EDOS to the accuracy 
necessary for this demonstration. Note that more rapid 
convergence can sometimes be achieved by generating k- 
point grids centred about a k-point other than T. We did 
not test such grids in this work. 

The results of these convergence tests are displayed 
in Fig. [7| Close examination of the enlarged EDOS at 
t = ps demonstrates that the occupied EDOS (below 
zero in energy) remains essentially unchanged with k- 
point sampling. The only noticeable effect is some reduc- 
tion in the sharpness of features when using the T-point 
with small numerical broadening (0.05 eV in this case). 
However, we see clearly that the unoccupied EDOS is 
greatly modified as the k-point density is increased. The 
gap that exists between LUMO and LUMO+1 (from 2 eV 
to 3.5 eV) under the T-point approximation is filled com- 
pletely at higher k-point densities. Furthermore, the 
qualitative form of the EDOS beyond the conduction 
band minimum is completley different. Comparison with 
the other uncorrelated molecular configurations reveals 
the same behavior. In fact, as a function of time, the 
converged EDOS shows more variation below the Fermi 
energy than it does in the unoccupied subspace. These 
variations are likely due to particular relative orientations 
of water molecules or making and breaking of hydrogen 
bonds, but will not be investigated here. 



VIII. CONVERGENCE OF THE EDOS WITH 
SYSTEM SIZE 

For an electronic structure calculation of a periodic 
system, accurate Brillouin zone integration may be 
achieved either by (i) increasing the k-point sampling 
of the Brillouin zone of the primitive unit cell, or (ii) 
by using larger supercells comprising repeated unit cells 
and a minimal k-point sampling. The latter approach is 
clearly more expensive and where possible we would pre- 
fer the former option. However, for disordered systems, 
the equivalence of these two approaches no longer holds, 
since the system is no longer periodic. Therefore, it is, 
in principle, more accurate within calculations performed 
under periodic boundary conditions to approach the limit 
of the bulk, disordered phase by increasing the supercell 
size and using the T-point approximation. This is also 
a requirement for the analysis of long-range molecular 
structure in the disordered phase. 

However, experience from tight-binding calculations of 
a variety of systems, indicates that the electronic struc- 
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FIG. 7: Convergence of the EDOS, computed within DFT/PBE, with respect to k-point sampling for 10 uncorrelated config- 
urations of 32 water molecules sampled every 20ps from a 200 ps TIP4P, NVT molecular dynamics trajectory at T=300 K, 
p = 0.996g/cm 3 . The EDOS are broadened using Lorentzians with a full width at half maximum of 0.05 eV. 
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FIG. 8: Comparison of EDOS for liquid water, computed within DFT/PBE, from 256 molecule supercells, sampled every 20ps 
from a 140 ps TIP4P, NVT molecular dynamics trajectory at T=300 K, p = 0.996g/cm 3 , with a time averaged EDOS from an 
uncorrelated TIP4P trajectory using a 32 molecule supercell, computed using 8 k-points in the first Brillouin zone. The EDOS 
are broadened using Lorentzians with a full width at half maximum of 0.05 eV. 
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ture of homogeneous systems is, perhaps, influenced most 
by short-range atomic structure, up to second-nearest- 
neighbor interactions. Therefore, we used small, 32 
molecule supercells and increased k-point sampling to 
approximate the electronic structure of larger disordered 
structures. Such comparisons are consistent when the k- 
point density is the samei^&Si This guarantees that the 
electronic degrees of freedom encompass the same vol- 
ume in real-space, the difference being that the system 
with larger k-point density possesses more structural or- 
der, and so is a poorer approximation to the disorderd 
liquid. 

In Fig. |H1 we present the EDOS computed using 
DFT/PBE for large, 256 water molecule supercells. 
These eight uncorrelated configurations of liquid water 
are taken from another independent TIP4P trajectory, 
prepared as outlined in Section[V] These EDOS are com- 
puted using the T-point approximation and compared 
with an EDOS representative of a 32 molecule liquid wa- 
ter supercell using 8 k-points. This EDOS for the smaller 
system is actually the average of the 8-kpoint EDOS over 
the 180 ps presented in Fig. [7| We see that the essential 
features of the EDOS of the larger system are well repro- 
duced by the smaller system with an equivalent k-point 
density. In fact the agreement is remarkable considering 
that there is no correlation between these structures and 
that the data spans 140 ps. 

The implications of this agreement are that smaller 
systems with increased k-point density can accurately 
represent the electronic structure of larger liquid water 
supercells. Furthermore, since we require at least 64 k- 
points to converge the EDOS of the 32 molecule configu- 
rations (Fig. UJ, it is clear that the 256 molecule, T-point 
calculations are still not converged. The computational 
expense of generating DFT/GGA electronic structures 
for such large systems prohibits the comparison with su- 
percells containing 2048 water molecules. Fortunately, it 
seems that such exhorbitant calculations are unnecessary 
for accurate prediction of (at least) the DFT electronic 
stucture of liquid water. 

We note also that the qualitative features of the 
Kohn-Sham eigenstate at the conduction band minimum 
(Fig. |5J) of these very large, 256 molecule supercells are 
the same as those of the LUMO computed in the smaller 
T-point calculations and in previous work. 



IX. CONSEQUENCES FOR SOLVATION AND 
SPECTROSCOPY 

At least for the electronic ground state of molecular 
systems, our calculations show that T-point sampling is 
sufficient for convergence of the DFT electronic structure. 
However, care should be taken to verify that the disper- 
sion in such systems is in fact minimal in order to justify 
this approximation. In particular, simulations involving 
phase transitions from disordered to ordered phases or 
from dilute to concentrated phases may exhibit different 




FIG. 9: An isosurface (gold) of the probability density of 
the Kohn-Sham eigenstate at the conduction band minimum 
of a 256 molecule supercell of liquid water extracted from a 
TIP4P trajectory. Water molecules indicated as red (oxygen) 
and white (hydrogen) ball-and-stick models. 



degrees of dispersion in the occupied EDOS, and the min- 
imum k-point sampling required for the more dispersive 
phase should be adopted. 

The large degree of dispersion in the lowest conduction 
band of liquid water (Fig. may have consequences in 
the simulation of a hydrated, excess electron. — In the 
limit of large system size, there should be a continuum of 
states at the conduction band minimum, allowing for the 
possibility of "intraband" transitions mediated by finite 
temperature. It is not clear what impact such transitions 
may have on the dynamics of such a system and the time 
scale for localization of the solvated electron. 

Of course, all of our work is limited by the accuracy 
of DFT, and, for excited electronic states, DFT has well- 
recognized limitations. This has been explicitly demon- 
strated for ice Ih in a recent publication^ However, 
many-body approaches, such as GW, 55 which allow for 
improvements in the description of the spectrum of ex- 
citations, rely on the use of Kohn-Sham eigenstates as a 
starting point. Therefore, issues concerning the super- 
cell size and k-point sampling are also relevant for these 
calculations. In particular, if the size of supercell can be 
reduced in favor of increased k-point sampling, this will 
be extremely advantageous for these computationally ex- 
pensive methodologies. 

DFT estimates of the optical absorption spectrum of 
liquid water will also be sensitive to the system size. We 
illustrate this in Fig. 1101 where we compare the joint den- 
sity of states (JDOS) of water for various k-point sam- 
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FIG. 10: Convergence of the joint density of states (JDOS) 
of liquid water (arbitrary units): for a 32 molecule supercell 
with 1 k-point (red); 8 k-points (blue); 64 k-points (black); 
and for a 256 molecule supercell with 1 k-point (green) with 
data only up to 10 eV excitations due to a limited number of 
unoccupied states in this calculation. A Gaussian broadening 
of 0.05 eV has been employed to expand the JDOS associated 
with each specific transition. 
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pling schemes and system sizes. The JDOS is a first es- 
timate of the optical absorption spectrum, ignoring the 
role of symmetry in electronic transitions from valence to 
conduction band. We compute the JDOS by considering 
only interband transitions from the valence to the con- 
duction band. We notice that the JDOS of liquid water 
is less sensitive to k-point sampling than the underly- 
ing EDOS (at least for the broadening scheme we have 
adopted in Fig. I10|) . This is to be expected given that 
the JDOS is a convolution of the valence and conduc- 
tion band EDOS. Furthermore, the same transferability 
to larger systems is apparent when we compare the JDOS 
of a 256 molecule liquid water supercell with that of an 
uncorrelated, 32 molecule supercell using 8 k-points. We 
note that both of these JDOS are qualitatively different 
from that computed using the smaller supercell within 
the T-point approximation. The existence of a peak at 
the absorption onset is one of the expected consequences 
of poor Brillouin zone sampling. 

DFT investigations into the optical properties of 
molecules and ions in aqueous solution also require par- 
ticular attention with respect to accurate representations 
of the electronic structure of water. Let us ignore, for 
the moment, the strong possibility of different system- 
atic band gap errors for the solute and solvent within 
DFT, as has been demonstrated using hybrid exchange 
correlation functional 56 If there exists a hybridization 
of solute and solvent states, which would modify the op- 
tical properties, this may be inaccurately estimated using 
DFT under the T-point approximation, particularly if the 
relevant, optically active solute state mixes with the bot- 
tom of the conduction band of water. ^^^^ Under 
the T-point approximation this state can hybridize with 
only one water state, whereas our analysis indicates the 
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FIG. 11: The "band structure" of a 32 molecule liquid water 
supercell in the presence of a core (x-ray) excitation from the 
Is orbital of one particular oxygen atom modelled using a 
modified pseudopotential (see text). 

presence of a continuum of states in this energy range. 

Finally, we demonstate (Fig. Ill|l the degree of disper- 
sion in the Kohn-Sham eigenstates of a 32 molecule su- 
percell of liquid water in the presence of an x-ray exci- 
tation. We model this excitation using a modified pseu- 
dopotential which includes a core hole in the Is level for 
one particular oxygen atom in the system. The system 
is effectively ionized and the impact of this perturba- 
tion is apparent when the band structure is compared 
with the ground state (Fig. 01. Localized states are re- 
alized, shifted from their respective bands. The conduc- 
tion band now exhibits states of varying localization, i.e. 
a subset exhibiting minimal dispersion and the remain- 
ing states retaining the typical dispersion associated with 
the ground state conduction band. Increased k-point 
sampling beyond the T-point is required to accurately 
describe all of the conduction band which will be incor- 
porated in spectral calculations for this system. 



X. COMPUTATIONAL EFFICIENCY 

The calculations presented here provide useful guide- 
lines, in general, for the estimation of the electronic prop- 
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erties of disordered molecular systems in the condensed 
phase. While these facts may be well known, we feel 
that it is useful to reiterate them here as they should be 
applied to liquid water. 

Firstly, it is important to know, for a given supercell 
size, what amount of k-point sampling is sufficient to de- 
scribe the occupied subspace of Kohn-Sham eigenstates. 
This can be easily checked by analysis of the band struc- 
ture of a representative molecular configuration, or, more 
consistently, by checking the convergence of the EDOS 
with respect to k-point sampling. Note that special care 
is required for the case of possible phase transitions, in 
small supercells, to states where dispersion is large. 

Secondly, as we have shown in this work, knowledge of 
the degree of dispersion in the occupied subspace is not 
necessarily transferable to the unoccupied subspace. It 
may be that a more dense k-point sampling is required 
for the unoccupied bands. However, given that these 
states do not influence the electronic charge density of 
the system, they may be generated non-self-consistcntly. 
One may quickly converge the electronic charge density 
of the system by considering only the occupied states and 
use a sparse k-point sampling, e.g., just the L-point for a 
32 molecule supercell of water. This charge density may 
be used to generate the common effective potential in a 
large set of Kohn-Sham equations which are only coupled 
if they correspond to the same k-point in the Brillouin 
zone. These equations may be solved using just one ma- 
trix diagonalization each, a process that is trivially par- 
allclizable with perfect linear scaling with respect to the 
number of k-points. 

Thirdly, we have shown that the EDOS of a rela- 
tively small supercell of water molecules computed with 
k-points is a very good approximation to the EDOS of a 
larger system, with an equivalent k-point sampling den- 
sity. In the case of water, this is a huge saving in com- 
putational cost, since increasing the system size eightfold 
introduces an increase in computational time by a factor 
of 512 for a typical planewave pseudopotential calcula- 
tion, which scales at worst as 0(N 3 ), where N is the 
number of electrons in the system. On the other hand, 
using 8 k-points for the original supercell introduces only 
an eightfold increase in computational cost, as we have 
shown. Furthermore, we have also shown that such large 
calculations: 256 molecules with T-point sampling; are 
still not converged in the EDOS. So, if we wish to ap- 
proach a converged result, increased k-point sampling 



may be the only resort. This necessity does not impose a 
limit on accuracy in electronic structure calculations for 
liquid water. 



XI. CONCLUSIONS 

In this work, we have provided a clear insight into 
the electronic structure of liquid water within the con- 
text of density functional theory - in particular using 
the PBE, gradient corrected, exchange correlation func- 
tional. We have shown that an accurate representation 
of the electronic structure is provided by relatively small, 
32 molecule supercells of the liquid. By inspection of the 
convergence of the electronic density of states (EDOS) 
with respect to k-point sampling of the first Brillouin 
zone, we have verified that the T-point approximation is 
adequate for accurate estimation of the occupied EDOS 
and consequently all ground state electronic properties. 
However, we find that the T-point approximation is inad- 
equate in providing an accurate description of the unoc- 
cupied EDOS of liquid water, which required a 4 x 4 x 4 
k-point mesh for convergence. Comparison with larger, 
uncorrelated, 256 molecule supercells of liquid water in- 
dicates that dense k-point sampling for small supercells 
provides an accurate estimation of the electronic struc- 
ture of larger supercells which more closely approximate 
the structural disorder of the liquid. This reveals the pos- 
sibility of marked savings in computational cost when ex- 
amining the electronic structure of molecular liquids such 
as water and, in particular, when attempting to use DFT 
to estimate the spectroscopic properties of such systems. 
Work is in progress to calculate both optical and x-ray 
absorption spectra for liquid water based on the analysis 
presented here. 
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